McmcHermes.jl

A documentation for the McmcHermes package.

McmcHermes provides a simple but efficient way to generate Markov Chain Monte-Carlo chains in order to sampling a probability density distribution.

Overview

The major functions in this module are:

one_mcmc: run a MCMC chain with one walker.

run_mcmc: run multiple chains with a specific number of walkers.

get_flat_chain: get the stored chain of MCMC samples.

get_gelman_rubin: get the Gelman Rubin convergence diagnostic of the chains.

Pkg Registry

using Pkg
Pkg.add("McmcHermes")

Example

First, let's generate some data:

using Distributions, Plots, LaTeXStrings, DataFrames, ProgressMeter

mu, sigma = 10, 2
l_b, u_b = 0, 20
d = Truncated(Normal(mu, sigma), l_b, u_b)
N = 1000
data = rand(d, N)

histogram(data, legend=false, size=(300,300), xlabel="data")

In order to sample the posterior probability distribution, it is necessary to define the likelihood, prior and logarithm of the posterior probability.

function log_likelihood(X::Vector, parameters::Vector)
    mu, sigma = parameters[1], parameters[2]
    y = 1 ./ (sqrt(2 * pi) .* sigma) .* exp.( -0.5 * ((X .- mu)./sigma).^2 )
    return sum(log.(y))
end

function log_prior(parameters::Vector)
    mu, sigma = parameters[1], parameters[2]
    if 5.0 < mu < 15.0 && 0.0 < sigma < 4.0
        return 0.0
    end
    return -Inf
end

function log_probability(X::Vector, parameters::Vector)
    lp = log_prior(parameters)
    if !isfinite(lp)
        return -Inf
    end
    return lp + log_likelihood(X, parameters)
end;

log_probability (generic function with 1 method)

Call the McmcHermes package and define the number of walkers, iterations, dimension of the parameter space and the initial guess.

using McmcHermes

mu, sigma = 10, 2
initparams = Vector{Float64}([mu, sigma])

n_iter, n_walkers = 1000, 200
n_dim, a = 2, 0.01

chain_tests = run_mcmc(log_probability, data, initparams, n_iter, n_walkers, n_dim, a=a);
size(chain_tests)

(1000, 200, 2)

Gelman-Rubin's diagnosis can be obtained from the chains.

println("Gelman Rubin Diagnostic: ", get_gelman_rubin(chain_tests))

Gelman Rubin Diagnostic: 1.0147180248917393

Finally, plot the chains.

labels = Tuple([L"\mu", L"\sigma"])
x = 1:size(chain_tests)[1]
p = []
for ind in 1:n_dim
    push!(p, plot(x, [chain_tests[:,i,ind] for i in 1:size(chain_tests)[2]], legend=false,
    lc=:black, lw=1, ylabel=labels[ind], alpha = 0.1, xticks=false))
end

plot(p[1], p[2], layout = (2,1))
plot!(size=(600,200), xlims = (0, size(chain_tests)[1]))

flat_chains = get_flat_chain(chain_tests, burn_in=100, thin=10)

flat = DataFrame(flat_chains, :auto)
colnames = ["mu", "sigma"]
flat = rename!(flat, Symbol.(colnames))

using PairPlots, CairoMakie
pairplot(flat)

• McmcHermes.jl
• • Overview
  • Pkg Registry
  • Example

Documentation for McmcHermes.